///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file input.h
///
/// \brief   header file for the class that holds user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_INPUT_HEADER
#define PCPS_INPUT_HEADER

#include<src/pcps.h>

namespace pcps {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   A class to hold user input options.
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  class Input {

    public:

      Input(int argc, char * argv[]);

      // functions to return the raw user input
      const std::set<std::string> & corr_sources() const { return _corr_sources; } ///< \sa _corr_sources
      bool verbose_print() const { return _verbose_print; } ///< \sa _verbose_print
      bool do_perturb() const { return _do_perturb; } ///< \sa _do_perturb
      bool properties() const { return _properties; } ///< \sa _properties
      bool fix_classical() const { return _fix_classical; } ///< \sa _fix_classical
      bool join_trotter_clusters() const { return _join_trotter_clusters; } ///< \sa _join_trotter_clusters
      bool random_correlators() const { return _random_correlators; } ///< \sa _random_correlators
      bool optimize_corrs() const { return _optimize_corrs; } ///< \sa _optimize_corrs
      bool optimize_agp_weights() const { return _optimize_agp_weights; } ///< \sa _optimize_agp_weights
      bool optimize_orb_coeffs() const { return _optimize_orb_coeffs; } ///< \sa _optimize_orb_coeffs
      bool optimize_pairing_matrix() const { return _optimize_pairing_matrix; } ///< \sa _optimize_pairing_matrix
      bool symmetric_pairing() const { return _symmetric_pairing; } ///< \sa _symmetric_pairing
      bool spin_to_fermion_dump() const { return _spin_to_fermion_dump; } ///< \sa _spin_to_fermion_dump
      bool project_sz() const { return _project_sz; } ///< \sa _project_sz
      bool precondition_sr() const { return _precondition_sr; } ///< \sa _precondition_sr
      bool sr_use_disk() const { return _sr_use_disk; } ///< \sa _sr_use_disk
      bool slightly_random_orbs() const { return _slightly_random_orbs; } ///< \sa _slightly_random_orbs
      bool active_energy_only() const { return _active_energy_only; } ///< \sa _active_energy_only
      bool s_half_orth() const { return _s_half_orth; } ///< \sa _s_half_orth
      bool spin_pair_correlations() const { return _spin_pair_correlations; } ///< \sa _spin_pair_correlations
      bool ssss_correlations() const { return _ssss_correlations; } ///< \sa _ssss_correlations
      bool spin_triple_products() const { return _spin_triple_products; } ///< \sa _spin_triple_products
      bool spin_triangle_3sz() const { return _spin_triangle_3sz; } ///< \sa _spin_triangle_3sz
      const std::string & boundary() const { return _boundary; } ///< \sa _boundary
      const std::string & correlator_type() const { return _correlator_type; } ///< \sa _correlator_type
      const std::string & hamiltonian_type() const { return _hamiltonian_type; } ///< \sa _hamiltonian_type
      const std::string & lattice_type() const { return _lattice_type; } ///< \sa _lattice_type
      const std::string & parallel_type() const { return _parallel_type; } ///< \sa _parallel_type
      const std::string & particle_type() const { return _particle_type; } ///< \sa _particle_type
      const std::string & ref_type() const { return _ref_type; } ///< \sa _ref_type
      const std::string & scalar_type() const { return _scalar_type; } ///< \sa _scalar_type
      const std::string & workdir() const { return _workdir; } ///< \sa _workdir
      const std::string & optimization() const { return _optimization; } ///< \sa _optimization
      const std::string & move_type() const { return _move_type; } ///< \sa _move_type
      const std::string & orb_source() const { return _orb_source; } ///< \sa _orb_source
      const std::string & agp_source() const { return _agp_source; } ///< \sa _agp_source
      const std::string & pairing_source() const { return _pairing_source; } ///< \sa _pairing_source
      double energy_tolerance() const { return _energy_tolerance; } ///< \sa _energy_tolerance
      double hubbard_t() const { return _hubbard_t; } ///< \sa _hubbard_t
      double hubbard_u() const { return _hubbard_u; } ///< \sa _hubbard_u
      double heisenberg_b() const { return _heisenberg_b; } ///< \sa _heisenberg_b
      double particle_spin() const { return _particle_spin; } ///< \sa _particle_spin
      double z_spin() const { return _z_spin; } ///< \sa _z_spin
      double damping() const { return _damping; } ///< \sa _damping
      double step_size() const { return _step_size; } ///< \sa _step_size
      double max_grad_norm() const { return _max_grad_norm; } ///< \sa _max_grad_norm
      double overlap_shift() const { return _overlap_shift; } ///< \sa _overlap_shift
      double cg_thresh() const { return _cg_thresh; } ///< \sa _cg_thresh
      double fci_shift() const { return _fci_shift; } ///< \sa _fci_shift
      double trial_energy() const { return _trial_energy; } ///< \sa _trial_energy
      double pmc_tau() const { return _pmc_tau; } ///< \sa _pmc_tau
      double tei_screening() const { return _tei_screening; } ///< \sa _tei_screening
      double init_config_thresh() const { return _init_config_thresh; } ///< \sa _init_config_thresh
      double uhf_desegragation() const { return _uhf_desegragation; } ///< \sa _uhf_desegragation
      double uhf_bias_decay() const { return _uhf_bias_decay; } ///< \sa _uhf_bias_decay
      int correlator_size() const { return _correlator_size; } ///< \sa _correlator_size
      int lattice_length_x() const { return _lattice_length_x; } ///< \sa _lattice_length_x
      int lattice_length_y() const { return _lattice_length_y; } ///< \sa _lattice_length_y
      int lattice_length_z() const { return _lattice_length_z; } ///< \sa _lattice_length_z
      int n_cor_orb() const { return _n_cor_orb; } ///< \sa _n_cor_orb
      int n_act_orb() const { return _n_act_orb; } ///< \sa _n_act_orb
      int max_hf_iter() const { return _max_hf_iter; } ///< \sa _max_hf_iter
      int max_cg_iter() const { return _max_cg_iter; } ///< \sa _max_cg_iter
      int max_cg_macro_iter() const { return _max_cg_macro_iter; } ///< \sa _max_cg_macro_iter
      int max_macro_loop() const { return _max_macro_loop; } ///< \sa _max_macro_loop
      int max_nr_iter() const { return _max_nr_iter; } ///< \sa _max_nr_iter
      int max_ac_length() const { return _max_ac_length; } ///< \sa _max_ac_length
      int nparticles() const { return _nparticles; } ///< \sa _nparticles
      int bfgs_length() const { return _bfgs_length; } ///< \sa _bfgs_length
      int dump_freq() const { return _dump_freq; } ///< \sa _dump_freq
      int print_freq() const { return _print_freq; } ///< \sa _print_freq
      int arena_block_size() const { return _arena_block_size; } ///< \sa _arena_block_size
      int arena_buffer_size() const { return _arena_buffer_size; } ///< \sa _arena_buffer_size
      int nthreads() const { return _nthreads; } ///< \sa _nthreads
      int global_seed() const { return _global_seed; } ///< \sa _global_seed
      int walkers_per_process() const { return _walkers_per_process; } ///< \sa _walkers_per_process
      int pop_control_freq() const { return _pop_control_freq; } ///< \sa _pop_control_freq
      int sr_disk_block_size() const { return _sr_disk_block_size; } ///< \sa _sr_disk_block_size
      long int sample_length() const { return _sample_length; } ///< \sa _sample_length
      long int burn_in_length() const { return _burn_in_length; } ///< \sa _burn_in_length
      long int properties_sample_length() const { return _properties_sample_length; } ///< \sa _properties_sample_length
      long int begin_true_projector() const { return _begin_true_projector; } ///< \sa _begin_true_projector
      const std::vector<int> & init_config() const { return _init_config; } ///< \sa _init_config
      const std::vector<double> & uhf_bias() const { return _uhf_bias; } ///< \sa _uhf_bias

      // functions to return data derived from user input
      int nsites() const { return _nsites; } ///< \sa _nsites
      int configs_per_site() const { return _configs_per_site; } ///< \sa _configs_per_site
      int n_agp_weights() const { return _n_agp_weights; } ///< \sa _n_agp_weights
      int n_rest_orb_coeff() const { return _n_rest_orb_coeff; } ///< \sa _n_rest_orb_coeff
      int n_pairing_variables() const { return _n_pairing_variables; } ///< \sa _n_pairing_variables

    private:

      // raw user input
      std::set<std::string> _corr_sources; ///< sources from which to build the initial correlators \sa pcps::keyword
      bool _verbose_print; ///< whether to print a lot or a little during execution
      bool _do_perturb; ///< whether to use perturbative long range pair correlators in the projected method
      bool _properties; ///< whether to compute additional wavefunction properties, beyond just the energy
      bool _fix_classical; ///< whether to fix the correlators generated by the "classical" correlator source
      /// \brief   whether to join together the clusters corresponding to different trotter
      ///          operators when computing the wavefunction update in time evolution
      bool _join_trotter_clusters;
      bool _random_correlators; ///< whether to randomize the non-fixed, non-ti correlators after they are read from their sources
      bool _optimize_corrs; ///< whether correlator variables should be optimized
      bool _optimize_agp_weights; ///< whether agp weight variables should be optimized
      bool _optimize_orb_coeffs; ///< whether orbital coefficient variables should be optimized
      bool _optimize_pairing_matrix; ///< whether pairing matrix variables should be optimized
      bool _symmetric_pairing; ///< whether to constrain the pairing matrix to be symmetric during optimization
      /// \brief   whether to, at the end of a spin calculation, convert the wavefunction to fermion form and dump it to a file
      bool _spin_to_fermion_dump;
      bool _project_sz; ///< whether to project the total spin of the wavefunction to a specific value
      bool _precondition_sr; ///< whether to use a diagonal preconditioner for the overlap matrix in stochastic reconfiguration
      bool _sr_use_disk; ///< whether to store the wavefunction derivatives on disk during the stochastic reconfiguration proceedure
      bool _slightly_random_orbs; ///< whether to add a little noise to the orbital coefficients after a Hartree-Fock calculation
      bool _active_energy_only; ///< whether to include nuclear and core energies when computing the local energy
      bool _s_half_orth; ///< whether to use S^(-1/2) orthogonalization when reading integrals from PSI3, otherwise reads the MO ints.
      bool _spin_pair_correlations; ///< when computing properties, whether to compute the spin dot products and Sz products on all site pairs
      bool _ssss_correlations; ///< when computing properties, whether to compute the four site Si.Sj Sk.Sl correlations
      bool _spin_triple_products; ///< when computing properties, whether to compute the icosidodecahedron's triangles' spin triple products
      bool _spin_triangle_3sz; ///< when computing properties, whether to compute the icosidodecahedron's triangles' sz*sz*sz products
      double _energy_tolerance; ///< energy tolerance used by the projected optimization method
      double _hubbard_t; ///< the hopping parameter in the Hubbard Model ( \f$ H_1 = -t * \sum_{\langle i j \rangle} a_i^+ a_j \f$ )
      double _hubbard_u; ///< the onsite repulsion in the Hubbard Model ( \f$ H_2 = U * \sum_i n_{i \uparrow} n_{i \downarrow} \f$ )
      double _heisenberg_b; ///< external field strength in the Heisenberg model in units of J
      double _particle_spin; ///< the spin of a single particle (e.g. 0.5 for fermions)
      double _z_spin; ///< if projecting to a specific total Sz, the value of Sz to project to
      double _damping; ///< damping applied to the first few updates of projected optimizations (larger damping means slower change)
      double _step_size; ///< step size used in various optimization methods
      double _max_grad_norm; ///< in steepest descent, gradients larger than this cutoff will be ignored
      double _overlap_shift; ///< shift to apply to the diagonal of the overlap matrix in stochastic reconfiguration
      double _cg_thresh; ///< threshold for solving linear equations with the conjugate gradient method
      double _fci_shift; ///< shift applied to the Hamiltonian's diagonal in the FCI method
      double _trial_energy; ///< trial energy used in fixed phase projector Monte Carlo
      double _pmc_tau; ///< timestep used in fixed phase projector Monte Carlo
      double _tei_screening; ///< screening threshold for the two electron integrals
      double _init_config_thresh; ///< threshold used for rejecting initial configurations with small weights \sa pcps::good_init_config_value
      double _uhf_desegragation; ///< parameter used to re-mix alpha and beta orbitals after a uhf calculation for the pfaffian initial guess
      double _uhf_bias_decay; ///< multiplier that reduces the uhf bias in each iteration of the uhf scf procedure
      int _correlator_size; ///< size of correlators automatically generated based on the lattice
      int _lattice_length_x; ///< extent of the lattice along the x dimension
      int _lattice_length_y; ///< extent of the lattice along the y dimension
      int _lattice_length_z; ///< extent of the lattice along the z dimension
      int _n_cor_orb; ///< number of core orbitals
      int _n_act_orb; ///< number of active orbitals
      ///< \brief   Maximum number of HF iterations used to generate initial orbitals.
      ///           If > 0, HF will be used to optimize the initial orbitals after they have been read from their source.
      int _max_hf_iter;
      int _max_cg_iter; ///< maximum number of conjugate gradient iterations allowed per solver attempt (i.e. each restart)
      int _max_cg_macro_iter; ///< maximum number of restarts for conjugate gradient solver
      int _max_macro_loop; ///< maximum number iterations for various optimization methods
      ///< \brief   Maximum number of Newton-Raphson optimization steps to apply in the projected optimization method.
      ///           Once these are exhausted, the method will revert to the sweep algorithm for any remaining iterations.
      int _max_nr_iter;
      int _max_ac_length; ///< maximum autocorrelation length to compute
      int _nparticles; ///< number of particles
      int _bfgs_length; ///< number of gradient vectors retained in the L-BFGS method
      int _dump_freq; ///< how frequently to dump the current wavefunction to a binary file
      int _print_freq; ///< how frequently to print the status of the optimization
      int _arena_block_size; ///< size, in bytes, of the memory blocks allocated by the arena memory manager
      int _arena_buffer_size; ///< size, in bytes, of the buffer allocated around blocks by the arena memory manager
      int _nthreads; ///< number of threads to use (not all optimization methods support threading)
      int _global_seed; ///< Overall random seed.  If < 0, the current time will be used instead.
      int _walkers_per_process; ///< target number of random walkers on each process for fixed phase projector Monte Carlo
      int _pop_control_freq; ///< frequency at which birth/death and overall population control is applied
      int _sr_disk_block_size; ///< size, in megabytes, of the blocks used for storing wavefunction derivatives in stochastic reconfiguration
      long int _sample_length; ///< total number of samples to be taken, this will be divided evenly between the processes/threads
      long int _burn_in_length; ///< number of burn in steps to apply to the markov chain on EACH thread
      long int _properties_sample_length; ///< number of samples to use when evaluating wavefunction properties other than the energy
      long int _begin_true_projector; ///< sample at which to switch to the true projector (i.e. no fixed nodes) in projector MC
      std::string _boundary; ///< the boundary conditions for the lattice \sa pcps::keyword
      std::string _correlator_type; ///< the type of correlator to build on the lattice (i.e. line, square, ...) \sa pcps::keyword
      std::string _hamiltonian_type; ///< the type of Hamiltonian to use \sa pcps::keyword
      std::string _lattice_type; ///< the lattice type \sa pcps::keyword
      std::string _parallel_type; ///< the way work is distributed across processes in the projected method \sa pcps::keyword
      std::string _particle_type; ///< the type of particle (e.g. fermion, spin, ...) \sa pcps::keyword
      std::string _ref_type; ///< the reference wavefunction that will be operated on by the correlator operators \sa pcps::keyword
      std::string _scalar_type; ///< type of numbers to use, either real or complex \sa pcps::keyword
      std::string _workdir; ///< the directory in which to write any files (WARNING, this defaults to ./)
      std::string _optimization; ///< the optimization method or type of calculation to apply \sa pcps::keyword
      std::string _move_type; ///< the type of move to use when evolving the markov chain \sa pcps::keyword
      std::string _orb_source; ///< where the initial orbital coefficients should be gotten from \sa pcps::keyword
      std::string _agp_source; ///< where the initial agp weights should be gotten from \sa pcps::keyword
      std::string _pairing_source; ///< where the initial pairing matrix should be gotten from \sa pcps::keyword
      std::vector<int> _init_config; ///< used to force the Monte Carlo walkers to start in a specific configuration
      /// \brief   The alpha vs. beta bias to apply during unrestricted hartree fock calculations.
      ///          In the hubbard model, this should be a single number representing the bias between up and down spin
      ///          occupations between the A and B sublattices.
      ///          For other hamiltonians, this vector should contain individual biases for each lattice site.
      std::vector<double> _uhf_bias;

      // data derived from user input
      int _nsites;  ///< the number of sites in the lattice (for fermions, alpha and beta sites are counted separately)
      int _configs_per_site; ///< the number of configurations a single site can take on (2 for fermions, 3 for spin = 1, etc.)
      int _n_agp_weights; ///< the number of agp weights
      int _n_rest_orb_coeff; ///< the number of restricted orbital coefficients
      int _n_pairing_variables; ///< the number of variables in the pairing matrix

      Input(); ///< the default constructor is disabled
      Input(const Input &); ///< the copy constructor is disabled
      Input &operator=(const Input &); ///< the assignment operator is disabled

      void verify();

      void initialize_derived_data();

      void process_input_command(const std::vector<char> & cmd, const int line_num);

      void bcast_input_command(std::vector<char> & cmd, int line_num);

      template<class T> void set_member(const std::string & member_name,
                                        std::stringstream & cmd_ss,
                                        T & member);

      template<class T> void set_member(const std::string & member_name,
                                        std::stringstream & cmd_ss,
                                        std::vector<T> & member);

      template<class T> void set_member(const std::string & member_name,
                                        std::stringstream & cmd_ss,
                                        std::set<T> & member);

      void set_member(const std::string & member_name,
                      std::stringstream & cmd_ss,
                      bool & member);

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   Holds keywords used by the PCPS program.
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  namespace keyword {

    // declare using static so that they have internal linkage
    static const std::string Part_Spin                            = "spin";
    static const std::string Part_Spinless_Fermion                = "spinless_fermion";
    static const std::string Part_Fermion                         = "fermion";
    static const std::string Boundary_Open                        = "open";
    static const std::string Boundary_Periodic                    = "periodic";
    static const std::string Lattice_Cartesian                    = "cartesian";
    static const std::string Lattice_Triangular                   = "triangular";
    static const std::string Lattice_Icosidodecahedron            = "icosidodecahedron";
    static const std::string Ham_Hubbard                          = "hubbard";
    static const std::string Ham_Heisenberg                       = "heisenberg";
    static const std::string Ham_File                             = "file";
    static const std::string Ham_PSI3                             = "psi3";
    static const std::string Ham_Molpro                           = "molpro";
    static const std::string Ham_Text_1212                        = "text1212";
    static const std::string Ham_Text_1122                        = "text1122";
    static const std::string Ref_Uniform                          = "uniform";
    static const std::string Ref_RHF                              = "rhf";
    static const std::string Ref_AGP                              = "agp";
    static const std::string Ref_Pairing                          = "pairing";
    static const std::string Ref_Pfaffian                         = "pfaffian";
    static const std::string Opt_Energy_Evaluation                = "energy_evaluation";
    static const std::string Opt_Steepest_Descent                 = "steepest_descent";
    static const std::string Opt_Stochastic_Imaginary_Time        = "stochastic_imaginary_time";
    static const std::string Opt_Stochastic_Real_Time             = "stochastic_real_time";
    static const std::string Opt_BFGS                             = "bfgs";
    static const std::string Opt_Projected                        = "projected";
    static const std::string Opt_Projected_TI                     = "projected_ti";
    static const std::string Opt_Imaginary_Time                   = "imaginary_time";
    static const std::string Opt_Real_Time                        = "real_time";
    static const std::string Opt_Full_CI                          = "fci";
    static const std::string Opt_PMC                              = "projector_monte_carlo";
    static const std::string Opt_Dump_DMRG_Ints                   = "dump_dmrg_ints";
    static const std::string Opt_None                             = "none";
    static const std::string Parallel_Site                        = "site";
    static const std::string Parallel_Farm                        = "farm";
    static const std::string Correlator_Line                      = "line";
    static const std::string Correlator_Square                    = "square";
    static const std::string Correlator_Triangle                  = "triangle";
    static const std::string Move_Standard                        = "standard";
    static const std::string Move_Exchange                        = "exchange";
    static const std::string Move_Mix                             = "mix";
    static const std::string Move_DoubleAB                        = "double_ab";
    static const std::string Scalar_Real                          = "real";
    static const std::string Scalar_Complex                       = "complex";
    static const std::string Corr_Source_None                     = "none";
    static const std::string Corr_Source_Binary                   = "binary";
    static const std::string Corr_Source_Text                     = "text";
    static const std::string Corr_Source_Classical                = "classical";
    static const std::string Corr_Source_Lattice                  = "lattice";
    static const std::string Corr_Source_All_Pair                 = "all_pairs";
    static const std::string Orb_Source_Binary                    = "binary";
    static const std::string Orb_Source_Text                      = "text";
    static const std::string Orb_Source_None                      = "none";
    static const std::string AGP_Source_Binary                    = "binary";
    static const std::string AGP_Source_Text                      = "text";
    static const std::string AGP_Source_RHF                       = "rhf";
    static const std::string AGP_Source_None                      = "none";
    static const std::string Pairing_Source_Binary                = "binary";
    static const std::string Pairing_Source_Text                  = "text";
    static const std::string Pairing_Source_Random                = "random";
    static const std::string Pairing_Source_AGP                   = "agp";
    static const std::string Pairing_Source_GHF                   = "ghf";
    static const std::string Pairing_Source_None                  = "none";

  }

}

#endif
